In this session we generalise the code from a growth model to a model with births and deaths.
There are many refinements to our code that we could make, and indeed we will work through several of them during the course, but our current framework is already sufficiently sophisticated to handle some much more complex population dynamics. We are therefore ready to start to thinking about our models rather than the mechanics of what we are doing.
In the modelling lecture we discussed a simple population growth model of the form
\[N(t+1) = N(t) + \lambda \times N(t)\]
which we can state roughly in English as
“the next population size is the current size plus the new additions (lambda times the current size)”
and this can be rewritten in turn in R as:
next.count <- latest$count + lambda * latest$countThis use of greek letters is fine for mathematicians, but we use growth.rate for lambda in our code as we want the code to read more like English. We will also write this on multiple lines as we calculate the new additions separately and then add them in. Note that, for consistency, throughout these practicals our code will refer to the data frame with the current latest population counts in as latest as in Practical 1-4, and these are referred to in the equations such as those on this page as \(N(t)\). The updated values, referred to above as \(N(t+1)\), appear as next.something in the functions if they have a name at all, and updated.something in the main script just to be distinct.
We now want to update the model to include both births and deaths, i.e.
\[N(t+1) = N(t) + b \times N(t) – d \times N(t)\] which can be similarly translated as
“the next population size is the current size plus the new births (b times the current size) minus the new deaths (d times the current size)”
where \(b\) represents the birth rate and \(d\) represents the death rate.
In this practical, we will reuse the RPiR package while adapting 0104-run-growth.R and 0104-step-growth.R appropriately. For convenience, we’ve generated 0105-step-birth-death.R (from 0104-step-growth.R) for you, with blanks for you to fill in:
#' ### Function: step_deterministic_birth_death()
#' Run one step of a simple deterministic exponential birth and death model.
#'
#' Arguments:
#' - latest -- the population count now
#' - birth.rate -- the birth rate
#' - death.rate -- the death rate
#'
#' Returns:
#' - the updated population count
#'
step_deterministic_birth_death <- function(latest, ___, ___) {
# Calculate changes to population
new.births <- ___ * latest$count
___ <- ___ * latest$count
# Calculate population at next timestep
next.count <- latest$count + new.additions - ___
# Create a data frame with the updated counts and return it
data.frame(count = next.count)
}## Error: <text>:12:52: unexpected input
## 11: #'
## 12: step_deterministic_birth_death <- function(latest, _
## ^Hint: Note that we give the names of the arguments to the function in the comments at the top.
#' ### Function: step_deterministic_birth_death()
#' Run one step of a simple deterministic exponential birth and death model.
#'
#' Arguments:
#' - latest -- the population count now
#' - birth.rate -- the birth rate
#' - death.rate -- the death rate
#'
#' Returns:
#' - the updated population count
#'
step_deterministic_birth_death <- function(latest, birth.rate, death.rate) {
# Calculate changes to population
new.births <- birth.rate * latest$count
new.deaths <- death.rate * latest$count
# Calculate population at next timestep
next.count <- latest$count + new.births - new.deaths
# Create a data frame with the updated counts and return it
data.frame(count = next.count)
}grade_this_code()## function (check_env)
## {
## check_env[[".__correct"]] <- correct
## check_env[[".__incorrect"]] <- incorrect
## grade_this(fail_code_feedback = fail_code_feedback, expr = {
## .message <- code_feedback(allow_partial_matching = allow_partial_matching)
## if (is.null(.message)) {
## pass(get(".__correct"))
## }
## else {
## fail(get(".__incorrect"))
## }
## })(check_env)
## }
## <bytecode: 0x7f8ec10ebd98>
## <environment: 0x7f8eb9df2120>Likewise, we’ve generated 0105-run-birth-death.R (from 0104-run-growth.R), with added blanks for you to fill in. Again,source() has been commented out as it won’t work in the {learnr} environment. Instead, the necessary files have been preloaded for you. Try to fill in the blanks and compare the outputs as you vary the birth and death rates:
library(RPiR)## Error in library(RPiR): there is no package called 'RPiR'step_deterministic_birth_death <- function(latest, birth.rate, death.rate) {
# Calculate changes to population
new.births <- birth.rate * latest$count
new.deaths <- death.rate * latest$count
# Calculate population at next timestep
next.count <- latest$count + new.births - new.deaths
# Create a data frame with the updated counts and return it
data.frame(count = next.count)
}library(RPiR)
# Load the step_simple_growth() function into the global environment (my workspace)
# Note this next line would be needed if you were running this in RStudio
# source(___)
#' Set up the simulation parameters
#' --------------------------------
#' First we set up the parameters for the simulation.
# Set the growth rate
human.annual.growth <- 0.015
# Set the death rate
___ <- ___
# Starting population size
initial.count <- 7000000000
# And setting times
start.time <- 0
end.time <- 100
#' Run the simplest possible simulation
#' ------------------------------------
#' Then run it so that we can get the output we need
# Set up the population starting size (at the first timestep)
population.df <- data.frame(count = initial.count)
# The timesteps that the simulation will run through
timesteps <- seq(from = start.time + 1, to = end.time)
# Now we loop through the time itself (starting at the second timestep)
for (new.time in timesteps) {
# Calculate population at next timestep
updated.population <- ___(latest = tail(population.df, 1),
___ = human.annual.growth,
___ = human.annual.death)
# Add new element onto end of population vector
population.df <- rbind(population.df, updated.population)
}
#' Plot the results
#' ----------------
#' And finally we output the results.
# Now we can plot the timesteps against the population vector
population.df$time <- c(start.time, timesteps)
plot_populations(population.df)
#' Set up a second population
#' --------------------------
___
#' Run the simulation
#' ------------------
___
#' Plot the results on the same plot
#' ---------------------------------
___## Error: <text>:14:1: unexpected input
## 13: # Set the death rate
## 14: _
## ^You can plot a second (or third, etc.) result on top of a first with commands like:
plot_populations(population.b, new.graph = FALSE, col = "red")However, if you want to do that, you should put the plot commands together, because the #' comment closes the current plot and you can’t then draw on top of it. Remember that col sets the colour of the line in the plot – use names “red”, “blue”, etc.